# This script will read in AWAP data, having already been extracted from Raw format by Lukes script
# It will then resize the files and change them to the appropriate datum.

# Establish directories

in.dir = '/home1/99/jc152199/'
setwd(in.dir)
out.dir = '/home1/99/'

# Load the gbm library to perform Boosted Regression Tree Analysis
necessary=c('SDMTools','sp','rgdal','raster')
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)
#load the libraries
for (lib in necessary) library(lib,character.only=T)

# Read in base ASCII file that is datum WGS1984-UTM, 250 m res, and sized to WTplusBuffer

#dem.asc = read.asc('dem.asc')

# Establish row/column positions for 250 m grid

coastdist.asc = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/coastdistwtplusbuffer250m.asc')
base.asc = coastdist.asc
#base.asc[,] = NA
base.pos = as.data.frame(which(is.finite(base.asc), arr.ind = T))
#west.pos = as.data.frame(which(is.finite(demwest.asc), arr.ind=T)) 

# Add east and north to this data frame of row/column positions

base.pos$east = getXYcoords(base.asc)$x[base.pos$row]
base.pos$north = getXYcoords(base.asc)$y[base.pos$col]

# Convert east-north to lat-long
	### EPSG 32755 WGS1984 Lat Long
	### EGSG 4326 WGS1984 UTM
	### EPSG 4916 GRS1980 Spheroid
	### ITFR1994 4336
	### EPSG 4019 unspecified datum
	
	tout = as.data.frame(spTransform(SpatialPoints(cbind(base.pos[3:4]),proj4string=CRS("+init=epsg:32755")), CRS("+init=epsg:4326")))
	tout2 = as.data.frame(spTransform(SpatialPoints(cbind(base.pos[3:4]),proj4string=CRS("+init=epsg:4326")), CRS("+init=epsg:32755")))
	
	
	base.pos$east = tout[,1]
	base.pos$north = tout[,2]
	
	maxlat = range(tout[,2])
	maxlong=range(tout[,1])
	
t.pos = data.frame(base.pos$col, base.pos$row, extract.data(cbind(base.pos$east,base.pos$north),dem.asc))

base.asc[cbind(base.pos$row,base.pos$col)]=t.pos[,3]

write.asc(x=base.asc, file=paste(out.dir,'dem250_ARCformat2.asc',sep=""))

#### Get north, south, east, and west from dem.asc

north = max(base.pos$north+125)
south = min(base.pos$north-125)
east = max(base.pos$east+125)
west = min(base.pos$east-125)
